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1 Introduction 



Consider an initial value problem 



y'{x) = f{x,y{x)), 1/(0) =1/0 



(1) 



To approximate the value y{t) many numerical algorithms such Runge-Kutta methods make 
computations for a set of points chosen on the interval [0,t]. Usually the chosen points form 
an arithmetic series with stepsize h. As h decreases, the algorithms yield more precise results. 
However, in practice it is often needed to find the approximation of y{t) with a given precision 
and then the problem of choosing an adequate value for h arises. A naive approach when h 
is divided by 2 if the iteration of the algorithm produces an error larger than allowed usually 
yields satisfactory results. With this approach, all information from the previous iteration 
is discarded. It is possible to use that information to predict the value of h which produces 
an error of an allowed magnitude. Using an adaptive iteration for h it is possible to make 
the use of the algorithm more computationally effective. For the sake of simplicity, we will 
show how to apply the adaptive stepsize technique to Euler integration, even though similar 
methods exist for Runge-Kutta and other numerical algorithms. We will start with deriving 
an adaptive recurrence relation for h from estimating the error Euler integration produces, 
and them compare the computational efficiency of the adaptive and naive approaches. 



Consider an initial value problem of the same form as (1). Suppose we arc interested in 
approximating y at a point x = t. For this purpose, split the interval [0, t] into n intervals of 
length h = -. Due to the intermediate value theorem for every x we have 



2 Standard Euler Integration 



y{x + /i) = y{x) + hf{x, y{x)) + ch? 



(2) 



where for some ^ e {x,x + h) we have 



(3) 
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This observation leads to Euler integration, a simple numerical method of solving ordinary 
differential equations. Let tk — kh ion: k — 0,1, ... ,n. Then define a recursive sequence as 
follows 

Vk+i = yk + hfitk, Vk), k>0 (4) 

Due to (2) the error introduced on each step of (4) is ch^ where c is proportional to y"{^). In 
practice, we most often need to approximate the value of y{t) with a given precision. If we use 
the standard Euler integration method with a fixed stepsize h, we cannot deduce anything 
about the error factor c. Therefore, we need to iterate the algorithm for larger values of n 
until the resulting ?/„ approximates y{t) with the desired precision. In other words we need 
to devise a sequence hi and compute the value of i/n^ using the relations i/q^ = yo and 

y^li-y'^^ + hif{h,y'i^) (5) 

for i = 1, 2, . . . until {yn^l^ — yn]\ becomes small enough. If the implementation of the algo- 
rithm uses fixed-precision representation of real numbers, which is often the case in practice, 
the complexity of each step of (5) is constant. Therefore, the complexity of computing yn\t) 
with an n— step integration is linear in n. Let rimin be the smallest number of steps necessary 
to achieve the desired precision. Then our goal is to reduce the total number of iterations 
of (5) and avoid values of n that are much larger than rimin- A natural way is to double 
the number of steps n on every iteration. In other words, hi is given by a simple recursion 
h,+i = |. This way only log2 n^m iterations of the algorithm are needed and all n that we 
use do not exceed 2n. 



3 Adaptive Stepsize Method 

It is easy to notice that every iteration of the Euler integration algorithm provides much in- 
formation that the simple approach when n is doubled on every step does not use. Moreover, 
the desired precision level is not used when calculating the stepsize on the next iteration. 
One solution to this problem which allows to optimize the iteration algorithm is to look at 
the error each step produces. We already found that each step of the iteration (5) produces 
an error of c/i^ where the rate of change of c is proportional to f"'{^). Let us define the error 
of each step 

4^ - vit) - yf (6) 

Then we are looking for |en^| < e to hold. If we assume that the error factor c in (2) is con- 
stant, it will enable us to devise a sequence hi that will lead to a more efficient approximation 
of y{t). If c is constant, then the accumulated error when calculating yn^ is 

^n- — '^ich^i — T chi — tchi — ahi (7) 

' hi 

where a is a constant. Consider two iterations with stepsizes hi and /i2- Recalling the 
definition of the error of each iteration, we have 

vit) - yS = = (8) 
vit) - yS - ^i? = (9) 
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Subtracting, we have 

z/i? - yS = - h,) (10) 

(2) _ (1) 

Let us find what is the condition for needed for the next iteration to produce an error 
within the desired range. 

^ I (3)1 I u I '''3 i/n2 yni 

^> 14/1 = ^31 = ^^^-3^ (12) 
{hi-h2)e 

^ |J2)_ m, ^13) 

This leads to a simple adaptive stepsize algorithm determined by the sequence hi which is 
defined as 

{hi-hi+i)e 

^ ^^^^ 

for some coefficient q < 1. 



4 Comparison 

Let us compare the efficiency of Euler integration using the naive doubling iteration and the 
adaptive stepsize method shown by applying them to the following initial value problem 

y'{x) = l-Ax + y{x), y{0) = 1 (15) 

The goal is to approximate the value y(l) with error less than lO""'. Wc will use the following 
simple implementation of Euler integration in Common Lisp using 8-byte floating point 
numbers. 



(defun euler ( f tn n yO) 

" Approximate ^y ( tn ) ^ wit h ^n — step ^ Euler ^method" 
(declare (type double — float x yO) 

(type integer n)) 
(let ((h (/ tn n)) (y yO ) ) 
(progn 

(loop for i from to (— n 1) 

do (incf y (* h (funcall f (* i h) y)))) 

y))) 



The exact solution of the initial value problem (15) is easy to find to be 



Then y(l) ^ 64.897. Let us choose the initial stepsize hi — 0.1 and find the values of yn} for 
the sequence hi defined by /ij+i = 2/ij first. 





34.411 


(10) 
t/c^i on 


^ 64.796 


(2) 


^ 45.588 


(11) 


^ 64.847 


(3) 


53.807 


(12) 


Ri 64.872 


(4) 

i/so 


^ 58 Ql fi 


(13) 
1/40960 


~ P,A 88 c; 


(5) 
1/160 


^ 61.786 


(14) 
1/81920 


64.891 


(6) 
1/320 


fti 63.310 


(15) 
1/163840 


64.894 


(7) 
1/640 


64.095 


(16) 
1/327680 


64.896 


(8) 
I/128O 


64.494 


(17) 
1/655360 


Ri 64.897 


(9) 
1/2560 


R! 64.695 







Therefore, 17 iterations are needed and the cost of computing the approximation is 

10m + 20m + • • • + 655360m = 10(2^^ - l)m = 1310710m (17) 

where m is the computational cost of one step of (5) which is constant for our implementation. 
Now let us keep the initial stepsizes hi — 0.1, /i2 = 0.05 but use the adaptive recurrence (14) 
with the coefficient q = 0.9. The first two approximations to y{t) are the same. 

~ 34.411 y^J ^ 45.588 
Because the allowed error e = 0.001, we have 

^ (0.1 - 0.05) X 0.001 , fi 
/i3 = 0.9^ Ri 4.026 X 10~^ (18) 

145.588 - 34.4111 ^ ^ 

Then na = \-^~\ = 248378. As expected, the iteration yields error close to the allowed limit. 

^248378 ~ 64.896 (19) 



For the next iteration we have 



/i2 - ^3 



h4 = 0.9——^ Ri 2.331 X 10-^ (20) 

164.896 - 45.5881 ^ ^ 

Then = f^] = 429086. The fourth iteration yields an approximation with the desired 
precision. 

yi?9086 ~ 64.897 (21) 

Thus, when a simple adaptive stepsize method is used, the number of iterations needed is 
only 4 with the total cost of computation 

10m + 20m + 248378m + 429086m = 677494m (22) 

which is around half of the cost of computing the approximation with doubling the number 
of steps. 
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5 Conclusion 

We have seen how to successfully apply the adaptive stepsize methods to Euler integration 
making it more computationally effective. Similar but more advanced techniques can be 
applied to more efficient numerical methods such as Runge-Kutta to develop adaptive step- 
size algorithms such as Runge-Kutta-Fehlberg and Dor mand- Prince methods which arc used 
in practice. For example, Dormand- Prince method is used in one of the Matlab ordinary 
differential equation solvers. 
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